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We consider the critical system with a point defect and study the variation of thermodynamic 
quantities, which are the differences between those with and without the defect. Within renormal¬ 
ization group theory, we show generally that the critical exponent of the internal energy variation is 
the specific heat exponent of a pure system, and the critical exponent of the heat capacity variation 
is that for the temperature derivative of specific heat of a pure system. This conclusion is valid for 
the isotropic systems with a short-range interaction. As an example we solve the two dimensional 
Ising model with a point defect numerically. The variations of the free energy, internal energy and 
specific heat are calculated with bond propagation algorithm. At the critical point, the internal en¬ 
ergy variation diverges with the lattice size logarithmically and the heat capacity variation diverges 
with size linearly. Near the critical point, the internal energy variation behaves as In |t| and the heat 
capacity variation behaves as where t is the reduced temperature. 

PACS numbers: 75.10.Nr,02.70.-c, 05.50.-l-q, 75.10.Hk 


I. INTRODUCTION 

As we know there are many dramatic effects in the 
critical phenomena because of the divercence of the cor¬ 
relation length at the critical point [l|. For example, 
the specific heat has a logarithmic divergence for the two 
dimensional Ising model 0. Are there some dramatic 
effects near the critical point when we add a point de¬ 
fect into the system? The answer is yes. The effects are 
dramatic and surprising. 

Consider the two dimensional Ising model with a point 
defect. We focus on the variations, the difference between 
with and without the defect, of thermodynamic quanti¬ 
ties. Naively, one would expect the internal energy vari¬ 
ation SU (X uq, where ug is the energy density of pure 
system and is a constant at the critical point. However 
in our numerical calculation on a. N x N lattice, 
at the critical point the internal energy variation 
is proportional to \ogN. Moreover the heat ca¬ 
pacity variation is proportional to N rather than 
logA^. Hence the defect makes an unexpectedly large 
contribution. 

Using renormalization group (RG) theory, we studied 
the influence of the defect on the thermodynamic quan¬ 
tities. Due to the long range correlation near the critical 
point, the defect changes the energy density across the 
whole system. We find that the internal energy variation 
is an integration of energy correlation, which gives rise 
to the specific heat. Then the internal energy variation 
is proportional to the specific heat of the pure system. 
Consequently the heat capacity variation is proportional 
to the temperature derivative of specific heat of apure 
system. With operator product expansion (OPE) [3, 
and conformal field theory (CFT) [ 1 , @|, we show these 
conclusions for the two dimensional Ising model explic- 
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itly. Then we use the scaling theory and RG to extend 
to common critical systems. 

To test the above conclusions, we study the two di¬ 
mensional Ising model on a finite size lattice with a site 
defect using bond propagation algorithm (BPA) [MU- 
The variations induced by the defect of the free energy, 
internal energy and heat capacity are calculated. The 
BPA results verify the RG and scaling argument about 
the internal energy and heat capacity variation. 

In fact, the effect of defects on the critical phenomena 
has been studied 40 years ago [T2l - [T3| . However those 
studies focused on the average effect of many defects 
rather than a single defect. In addition, this study is not 
purely academic, but has potential practical use. The 
rapid development of Nano techniques makes it possible 
to study the critical phenomena on small size systems 
[IM3- In the production of the samples, the defect is 
unavoidable usually. Therefore it is of interest to study 
the effect of a defect in a finite size system. 

We arrange this paper as follows. In the section B, we 
present RG and scaling argument. In section IB, we re¬ 
port the BPA results on the two dimensional Ising model 
with a defect. Section IV is a concluding remark and 
acknowledgement. 


II. THE CRITICAL EXPONENTS OF THE 

INTERNAL ENERGY AND HEAT CAPACITY 
VARIATION 

We shall prove generally that for a spin system with 
short range interaction in any dimension: 

1. The critical exponent of the internal energy varia¬ 
tion due to a site defect is the same as that of the specific 
heat of a pure system. 

2. The critical exponent of the heat capacity varia¬ 
tion is the same as that of the temperature derivative of 
specific heat of a pure system. 

At first we show these conclusions for a two dimen- 
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To study the effect of the defect, we will calculate the 
variation of free energy defined by 

SF = Fo-F ( 6 ) 

where F = —InZ and Fq = — In^o; the internal energy 
variation 


SU = Uo-U (7) 

where U = dF/dp and 17o = dFo/dP; and the heat ca¬ 
pacity variation 


6C = Co-C ( 8 ) 

where C = d^F/dp"^ and Cq = d^Fo/dp^. 


FIG. 1. (a) 7 X 7 lattice with a defect in the centre. The 

spin 0 is labelled by the black spot at the center. Its four 
neighbors are labeled by 1,2, 3, 4. The bonds between spin 
0 and 1,2, 3, 4 are zero. All other bonds are equal to 1. (b) 
After applying Y — A transformation to the four spins in (a), 
where the center spin is spin 2, we get this lattice. (c)Using 
BPA, we get this lattice, (d) Applying the same procedure in 
three other directions, the original 7x7 lattice is transformed 
into a 6 X 6 lattice. 


sional Ising model explicitly. Then we give a general 
scaling argument for common critical systems. 

Consider the two dimensional Ising model on a square 
lattice with a site defect. The Hamiltonian is given by 

( 1 ) 

with Joi = Jo 2 = </o 3 = Jqa = 0 and all other bonds 
Jii! = 1. The summation is over the nearest neighbours. 
A 7 X 7 such lattice is shown in the Fig. 1. We will study 
this model with BPA in the next section, so we also show 
the schematic of the BPA. Its partition function is given 
by 


Z = Tre-^^ = ^ (2) 

{^i} 

The Hamiltonian for the pure system is given by 

Ho = - y] CTiCr^/. (3) 

and its partition function is given by 


A. The RG and CFT argument for two 
dimensional Ising model with a defect 


For the free energy variation, we have 


SF = In 


Tr{e-P^») 


= ln(e-/^-), 


( 9 ) 


Here we introduce the notation (• ■ • )o 
Tr[- ■ ■ . Using the equation 




X - Ya^aj, 


( 10 ) 


where X = cosh I3,Y = sinh /3, we can expand e 
into 

= Y[{X- Yaoa,) 

4 4 

= Ai - A2 y] cToCTj -b A3 y] cTiCTj 

i—l i<j 

4 

— H4 (TofTifJjCTfc + A5(Ti(J2<T3Cr4 

i<j<k 

( 11 ) 

where the coefficients are given by Ai = A"^,A 2 = 
X^Y, A 3 = X^Y^,A 4 = XY^, As = The free energy 
variation is given by the average of these local variables. 
It is singular since Fq is singular. We will show the singu¬ 
larity of the internal energy (derivative of SF) variation 
explicitly. 

For the internal energy variation we have 


Zo = Tre-^^° = (4) 

{<^i} 


Uo 


Tr{Hoe-^^o) _ , ^ 

rr(e-/5ffo) 


( 12 ) 


Define the variation of the Hamiltonian 


and 


4 

(5H = H - Ho = y] ao(T,. 


u = 


(e-P^^iHo+SH))^ 

(e-/5^^)o 


( 5 ) 


(13) 
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Then we get 


5U = 




-{H^) 


0 • 


(14) 


Substituting equation (HU into 5U, we get an expan¬ 
sion. Because is a higher order 

contribution as compared to (e“^^^iJo)g / j 

ignore it. The contribution to SU from the second term 
in Eq. m is given by 


4 

E E Q {^k^k')(j)- (lb) 

2—1 <k,k'> 


Here and below we neglect the denominator 
which is a finite quantity and can be replaced by its value 
at the critical point. The quantity in the above equation 
is proportional to the specific heat of the pure system. 
The second non zero term in the expansion oi SU related 
to the third term in Eq. (HU is given by 

4 

E E {{(JiCTjO'kO'k'') Q Q {^k^k^} {)) • ( 1 ^) 

i<j <k,k'> 

With OPE in RG we will show both terms in the above 
two equations are of the same order and proportional to 
the specific heat of a pure system. Generally^ product 
of two nearby basic operators is reducible as H, d, d, 

Oa(ri)Oa(r2) = y^Ca/3,-y(ri2)l?('^(Rl2) 

7 

r = ri-r2, R=(ri-fr2)/2 (17) 

where Cap^-y are OPE coefficients. The sum is over all 
the scaling operators. The exact meaning of the above 
equation is defined by dl 


{Oa{yi)Oa{v 2 ) • • • )o = (r 12 )(^^ (Rl 2 ) ■ ’ ' )o 

( 18 ) 

The point is that, in the limit when |ri — r 2 | is much 
less than the separation between ri,r 2 and all the other 
arguments in • • •, the coefficients Cct/ 3,7 are independent 
of what is in the dots. 

For two nearby spin operators in two dimensional Ising 
model, we know the exact expansion from the conformal 
field theory (GET) @ 

-1 1 3 

CToo-i--|roi| 4 -f-|roihe(Roi)-I- (19) 

where e is the energy density operator, and • • • are irrel¬ 
evant operators with larger scaling dimensions. The en¬ 
ergy operator e(r) is also called thermal operator. Note 
that e is the operator for energy density minus its criti¬ 
cal value, so (e(r))o = 0. Then aiaj with 1 < f < j < 4 
in Eq. (1161) can also expanded as linear combination of 
identity and energy density operators. Both Eq. m 


and (HU) give rise to energy correlators. For a energy 
correlator, 

(e(0)e(r))g ^ ^ (20) 

where is the correlation length for t > 0 and t < 0 
respectively, and t is the reduced temperature.. 

From Eq. (HU), we get their contribution of terms in 
Eq. (ITSl) . (ITBl) to the internal energy variation 

E E ^|i’o2|^|i'fcfe'lbe(R'02)e(Rfcfe'))o 

i—1 <k,k'> 

+ E E ^l’'bl^kfcfc'lbe(R-2j)e(R-fcfc'))o 

i<j <k,k'> 

1 3 3 3 

^4 E i^fcfc'ivEi^o^i^i+Ei^bi^) ^2 

<k,k'> i—1 i<j ^ 

( 21 ) 

Considering that i,j are the neighbors of the defect and 
k,k' are neighbors, the approximation |Rfcfe' — Roi| ~ 
|R/c/c' — Ryl ~ Tk is used, where is the separation be¬ 
tween the spin k and the defect. One may note that in the 
summation of Eq. (HU, (HU) there are terms {o'iaj(Tk'^k')Q 
with the positions of ak,Uk' are not far from the defect. 
For such terms, we can not get the simple energy corre¬ 
lators in the above equations. However the main contri¬ 
bution comes from those terms with au , <Jk' far from the 
defect. For these terms, we can apply OPE to expand 
CTiCTj and UkCfk' with Eq. (USD. 

The left two terms in Eq. (ED are products of four 
operators. For the fourth term in Eq. ED, we take the 
term (Tocricr 2 cr 3 as an example. It can be expanded as 

(ToUifJ2U'3 — C(7(T(T(T,/ 4“ Ccrcr(T(T,€C (R 0123 ) + • • • ■ (22) 

In this expansion, another relevant scaling operator tr 
does not exist because it has different symmetry with 
the product of four spin operators. The coefficients 
Caaaaj, Caaaa,€ are functions of ro,ri,r 2 ,r 3 . The po¬ 
sition R 0123 = (ro + ri + r 2 + r 3)/4 is a result by consid¬ 
erations of symmetry. It is obvious that IR 0123 I < 1- The 
explicit form of Co-o-o-o-,/, Co-o-o-o-.e can be obtained from the 
calculation of 6-point correlation, say (o'ocri(T 20 ’ 3 crfc(Tfe')o, 
where (ikjO'k' are nearest neighbori ng spins and far away 
from the defect, i.e. |rfc|, |rfe 2 | » 1 Q. However we don’t 
need the explicit form of these coefficients here. 

Substituting it into Eq. ED , we get a energy correlator 
(e(Roi23))e(Ri/'))o l/’'z^> where I, I' are over the whole 

system. All the four spins terms in Eq. ED can be dealt 
similarly. Therefore all the terms in SU are proportional 
to the integration of the energy correlator. 

To be more explicit, we see the scaling further. At the 
critical point (e(0)e(r))g oc 1/r^, then the summation 
over the system leads to 

SU ^ ^ — oc log N (23) 

<ii'> 





4 


if the system size is N. Near the critical point, 
(e(0)e(r)p oc If TV the summation over 

the system gives rise to 

6U ^ ^ -log 1^1 (24) 

<w> 

This is just the behaviour the specific heat of a pure 
system. Then we have 

SU oc c (25) 

As we can see the influence of the defect on the inter¬ 
nal energy is an integration of energy correlator, which 
gives rise the fluctuation of energy. The physical nature 
of specific heat is just the fluctuation of energy. Due to 
the long range correlation, the defect changes the energy 
density across the system. The critical exponent of the 
internal energy variation is the same as that of the spe¬ 
cific heat of the pure system. 

Since SC = it has 

6C^-^\t\-\ (26) 

We simply get the heat capacity variation. Using finite 
size scaling [^, at the critical point, we get 

SC oc N. (27) 


B. The general scaling argument for an isotropic 
system with a defect 


Consider a general model with short range interaction 
in d-dimension. 


H = Ja^aj (28) 

<ij> 

The basic operators, are the spins, for example, Ui = ±1 
in Ising model, ai = ((^ix,criy),af^+afy = 1 in XY model, 
etc. 

The defect will break the bonds which connect to the 
defect. Generally we have 

g-/35it ^i + SH + ^ (SH)^ + ... = Y baSa (29) 

a 

where Sa are the all possible products of a^ai, i = 
1,2, 3,4 are the four nearest neighbours of the defect, 
and ha are the coefficients. According to the OPE, these 
products of basic operators can be expanded as a linear 
combinations of scaling operators. That is to say Sa can 
be expanded as 


Sa = ^Cak(l>k{^a) (30) 

k 

where (j)k are the scaling operators and Cak are the OPE 
coefficients . The position can be determined accord¬ 
ing to OPE in Eq. m step by step. Because the position 


Ti of the basic operators ai at the defect and its nearest 
neighbours satisfy jr^l < 1 if we set the lattice constant 
be the unit length, it should have |ra| < 1. 

For the scaling operator (j)k, the scaling dimension is 
Xk = d — Uk, where d is the spatial dimension and yk is 
eigenvalue of Uk in the RG. At the critical point the corre¬ 
lation between two scaling operators ^fc(O),0z(r) decays 
as 

(</.fe(0)</.z(r))o^(31) 

in the limit of r —>■ oo. For a critical system, the most 
relevant operator are magnetic operator s(r) and thermal 
operators (or energy operators) e(r). Their eigenvalues 
are yh and yt respectively and yh,yt > 0. The other 
scaling operators’ eigenvalues are negative and hence ir¬ 
relevant. 

In the expansion Eq. O, the magnetic operator s(r) 
is prohibited because of symmetry. We may classify the 
scaling operators as being even or odd under the sym¬ 
metry (Ti —^ — CTi. In the expansion Eq. (sni), the scal¬ 
ing operators (j)k('ra) must be even. Besides the constant 
term, the first order term is the thermal operator e(ra). 
Then, for the internal energy variation defined in Eq. 
d, the leading term is a summation of energy correla¬ 
tor between the defect and other bonds. The correlator 
of thermal operators near the critical point behaves 

— 

(e(0)e(r))o ^ J.2(d-yt) ~ J.2(d-i/v) ■ (^2) 

where v = 1/yt is the correlation length exponent and 
Ttie correlators between thermal operator 
and other even scaling operators decays faster since the 
scaling dimension of other even operators are larger than 
thermal operator’s. 

Therefore at the critical point for a finite system with 
size N and A^ ^ 1 , the internal energy variation is given 

by 

^;5^^A + TV‘^-2A = A + TV“/‘^ (33) 

where a = 2 — dv is the specific heat exponent, A is 
a constant, a = 0 is the logarithmic case. If zz > 0, 
the constant term A is negligible, and SU is divergent 
as A^ —>■ oo. Otherwise A is the leading term and SU 
converges to A as A^ —>■ oo. For a finite system the specific 
heat at the critical point just scales as c(L) A + . 

For away from the critical point and N ^ we have 


SU 



rd lfj:y 
j.2d—2jv 


A'+ - A" + |tr' 


(34) 


Similarly if z/ > 0, the constant term A" is negligible, 
and SU is divergent as |f| —>■ 0. Otherwise A" is the 
leading term and SU converges to A" as |t| —^ 0. In this 
case SU has a cusp. As we know for a pure system the 
specific heat scales as c ~ A" -I- |t|““ near the critical 
point. Therefore we have SU ^ c. 
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Since 5C = it has 


5c itr“-^ 


For a finite system at the critical point 


27| 


5C - 


(35) 


(36) 


The above argument is quite general except we as¬ 
sumed that the system is isotropic, i.e. the correlator 
only depends on the separation between two points. Thus 
conclusions in Eq. (l34l) and (1351) are valid for common 
isotropic systems. 


III. NUMERICAL SOLUTION OF THE TWO 
DIMENSIONAL ISING MODEL WITH A 
DEFECT 


To test the above conclusion, we calculate the internal 
energy variation and specific heat variation of two di¬ 
mensional Ising model with a defect numerically. In the 
lattice shown in Fig. 1, the boundary condition is open, 
i.e., we have four edges and four corners. For this kind 
of geometry, BPA is very efficient [MU- The sketch of 
BPA for this problem is also shown in Fig.l. To keep the 
defect be the centre of lattice, we set the size be odd and 
the aspect ratios are also odd. In this way we can avoid 
the trouble stemming from the asymmetry. 

BPA is very accurate and can be carried on lattices 
with extremely large sizes, say 2000 x 2000. With it, we 
have verified the CFT predictions on the corner free en¬ 
ergy with free boundary condition [a [2^. The result on 
the central charge of Ising model agrees with CFT in the 
precision of 10“^° Q. We also recover the aspect ratio 
dependence of the corner free energy in CFT theory accu¬ 
rately The corner free energy with fixed boundary 

and mixed boundary condition has been studied by CFT 
recently [211-1^1 . With an extended BPA , we verified 
these CFT prediction in the accuracy 10“^® |25|. 


In the algorithms, the transformations preserve these 
quantities during every step. The accuracy is only limited 
by the machine’s accuracy. With these algorithms, we 
can calculate the free energy, internal energy and specific 
heat with the same accuracy. As discussed in reference 
0, the time of calculation is proportional to x M 
if the lattice size is M x L. Then the accumulation of 
round-off error is proportional to L^/M. Therefore the 
accuracy can reach 10“^® — 10“^® for the lattice with a 
very large size 1000 x 1000 if all the variables are assigned 
in the quadruple precision format, in which the machine 
accuracy is 10“®®. This has been shown in EH- 
We apply the BPA to calculate the free energy and in¬ 
ternal energy directly. For the heat capacity, we adapt a 
difference approximation. In fact, the BPA for the heat 
capacity on the usual lattice has been developed [HI- 
However, for the lattice with a defect, we meet an un¬ 
known problem in the BPA for the heat capacity. There¬ 
fore, we calculate the internal energy U{f5 + A/3) and 


TABLE I. The the fitted Aq in Eq. (I37II for the free energy 
variation for p = 1,3, 5, 7, 11. It has Smax < 10“^^. 


p 

^0 

1 

1.409210310117490418345(2) 

3 

1.409210310117490418347(1) 

5 

1.409210310117490418347(1) 

7 

1.409210310117490418346(2) 

11 

1.409210310117490418344(5) 


U(/3 - A/3) and take [f7(/3 -h A/3) - f7(/3 - A/3)]/(2A/3) 
as the approximate value of dU/djS. In the numerical 
calculation, we set A/3 = 10“®® and set all variables in 
quadruple precision format. Because the result of BPA 
for a lattice with 1000 x 1000 can reach 10“^®, such a 
small difference A/3 = 10“®® does not cause serious in¬ 
stability. Then we can get the specific heat variation in 
a very high accuracy. 


A. At the critical point 

At first, we calculate the variations at the critical point 
/3 = /3c = 5 ln(l -I- V^). To investigate the geometrical 
effect, we study the rectangle with size MxN with the as¬ 
pect ratios p = M/N = 1, 3, 5, 7,11. As mentioned above 
to keep the defect be the centre of lattice, we set the size 
be odd and the aspect ratios are odd. For p = 1,3, 
we calculate 65 data points with 31 < A < 1021. For 
p = 5,7, 11, the number of data points are 63, 53,45 re¬ 
spectively and the range of size 31 < A < 981,781,621 
respectively. We fit the variations within finite size scal¬ 
ing. 

We fit the free energy variation with 

krnax 

SF=J2 ( 37 ) 

k=0 

In our fitting it has kmax = 12. This indicates that our 
numerical calculation is very accurate. The first term Aq 
gives the defect’s contribution to the free energy in the 
thermodynamic limit A —>■ oo. We show the first term for 
p = 1,3, 5,7,11 in Tab. I. In the thermodynamic limit, 
the (bulk) free energy density for a pure Ising model 
(without the defect) is given by fbuik = 0.929695398 • • • 
at the critical point 0. The contribution of the defect 
to the free energy is about Aq/ fbuik ~ 1.5 times of the 
bulk free energy density. The other fitted parameters 
can be seen in appendix A. The high order corrections 
stem from the boundary and finite size effect. Only Aq is 
boundary independent. That is to say, if someone solve 
this problem with periodic boundary condition, the Aq 
should be the same and the higher order corrections may 
be different. 

The fitting method is the Levenberg-Marquardt 
method [2^ for nonlinear fit. To characterize the ac- 
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TABLE 11. The fitted Bio in Eq. (1391) the internal energy 
variation for p = 1, 3, 5, 7, 11. It has Smax < 10“^^. 


p 

Bio 

1 

1.7560009405073355(2) 

3 

1.75600094050733606(1) 

5 

1.75600094050733605(1) 

7 

1.7560009405073360(2) 

11 

1.756000940507336(5) 


curacy of our fittings, we define the maximal deviation 

Smax = Max\yi - (38) 

where yi is the numerical data and y/** is the value given 
by the fitting formula. We choose the maximum of the 
deviations from the data to the fitted ones to represent 
our fitting quality. For a given fitting formula , the devi¬ 
ations are given by the nonlinear fitting algorithm. We 
expand the the free energy, internal energy and specific 
heat to as high order as possible to make the Smax as 
small as possible. In every table of the fitting parame¬ 
ters, we give the maximum of deviation. For example, in 
Tab. m Smax = 10-22. 

As mentioned above, the results of BPA for F, Fo,U, Uq 
are in an accuracy of 10-2®—10-2®. The fitting accuracies 
of SF,SU are about IO -22 — IO- 2 ®. They are consistent 
because F, Fq, U, Uq are in order of 10® for 1000 x 1000 lat¬ 
tice, the first six digits of F, Fq are cancelled in SF and it 
is the same for U, Uq and SU. This is why the maximal de¬ 
viations in Table I-X are 10-22 — 10-2®. For heat capacity 
variation, we use SC ~ [5C/(,5-l-A/3) —(517(/3—A/3)]/(2A/3) 
and A/3 = IQ-^^, then the accuracy of SC should be in 
order of 10-^^ — since the deviation of SU is mul¬ 

tiplied by l/A/3. In the fittings of SC, the maximal de¬ 
viation is about 10-®2 as shown in Tables XI-XV. This 
is also consistent. 

The internal energy variation can be fitted by the for¬ 
mula 

BiklnN + Bok , , 

»= E jv-t (39) 

k=0 

In our fitting, kmax is set to be 8. The leading term di¬ 
verges logarithmically with the system size. This leading 
term is obtained from Eq. (USD- The coefficient of leading 
term Bm is shown in Tab. II. As we can see the leading 
term is geometry independent, i.e. Biq are the same for 
p = 1,3, 5, 7,11 in the error range < IQ-^®. The higher 
order terms can be seen in the appendix B. Biq should 
be independent of the boundary condition. Therefore, if 
someone solve this problem with periodic boundary con¬ 
dition or other boundary condition, the Biq should be the 
same and the higher order corrections may be different. 

We fit the heat capacity variation in the following way 

= + + (40) 

fc =0 


TABLE III. The fiited C-i and C 20 in Eq. (1401) for the 
specific heat variation for p = 1,3,5,7,11. It has Smax < 
10-^2 


P 

C-i 

C20 

1 

0.370131893446(5) 

0.598836(1) 

3 

0.473858336337(2) 

0.59883789(1) 

5 

0.475120770039(2) 

0.59883783(2) 

7 

0.475126998263(7) 

0.59883752(4) 

11 

0.475127020038(2) 

0.59884(1) 


In our fitting, kmax is set to be 6. The leading term C_i 
and C 20 are shown in Tab. III. The heat capacity vari¬ 
ation shows a more dramatic effect. It diverges linearly 
with the size! This leading divergent term is obtained 
from Eq. (EZl) in the last section. The coefficients of this 
linear term depend on the aspect ratio p, i.e. they are 
geometry dependent. The next leading term is C 20 In^ N. 
Contrary to the first leading term, this term C 20 is inde¬ 
pendent of the aspect ratio, at least in the error range. 
The other fitted parameters can be seen in the appendix 
C. This fitting formula, especially for the term ln2 N, is 
a bit strange. We have tried many different formulas to 
fit the data. Only with this formula, we can fit the data 
with the smallest deviation. 

One may note that the geometry in our considera¬ 
tion as shown in Fig. 1 has open boundary: open edges 
and sharp corners. These boundaries will produce edge 
terms and corner terms [1) [13) nil- However, in our sit¬ 
uation, the boundary effects are cancelled exactly since 
both F,U,C and Fq,Uq,Cq contain the same boundary 
terms. We have expanded F, Fq separately and found 
that their edge terms are the same and so are their cor¬ 
ner terms. It is the same for U,Uq and C,Cq. Hence in 
SF, SU, SC there are no boundary terms. 

From the fitting at the critical point, we get SU ^ In A 
and SC ~ N. This verifies the conclusions given in Eq. 
(ESI and (EZD. 


B. Temperature dependence of the variations 

To study the finite size scaling, we calculate the varia¬ 
tions for square lattices with size N = 101,201,401,801. 
Fig. 2 shows the variation of the free energy in the criti¬ 
cal regime for square lattice with these sizes . In order to 
see the ratio between SF and the bulk free energy den¬ 
sity fbuik, we also show the bulk free energy density 0 
with black line in Fig. 2a. Near the critical point, the 
ratio SF/fbuik is about 1.5 for N > 100. As we can see 
these results almost coincide except for very small t be¬ 
cause the lattice sizes are very large. Fig. 2b shows the 
convergence of the results of different size at very small 
t. 

It is should be emphasized that we use directly the 
reduced temperature t = {T — Tc)/Tc, where = 
2/ln(l -I- v^) = 2.2691853- • • is the infinite volume Tc- 
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FIG. 2. (a) The free energy variation in the critical regime, 

(b) The finite size effect on the free energy variations. 
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FIG. 3. (a) The internal energy variation in the critical 

regime, (b) The logarithmic dependence of the reduced tem¬ 
perature for t > 0. (c)The logarithmic dependence of the 
reduced temperature for t < 0. The legend in (a) and (c) is 
the same as that in (b). 


Because the lattice sizes N = 101,201,401,801 are large 
enough so that the finite-size critical temeprature is very 
well described by the infinite volume Tc. 

To verify Eq. (IMl) and (|26l) with the the finite size 
scaling, we also study the dependence of the variations 
on the temperature. Fig. 3 shows the internal energy 
variation for the lattices with size N = 101, 201,401, 801. 
Fig. 3a shows the logarithmic divergence. In order to 
show this divergence more clearly, we show the part for 
t > 0 and f < 0 in Fig. 3b and 3c respectively. The linear 
dependence of the variation on the In |t| can be seen. As 
the size increases, the linear region grows. Fitting the 


FIG. 4. (a) The heat capacity variation in the critical regime, 
(b) The inverse dependence of the reduced temperature for 
t > 0. (c)The inverse dependence of the reduced temperature 
for t < 0. The legend in (a) and (c) is the same as that in 
(b). 


linear parts for N = 801, we get 

6Uf=iB'ln\t\ (41) 

where B' = 1.888(4) for 0.003 < t < 0.1 and B' = 1.68(5) 
for 0.003 < |t| <0.1 and t < 0. This agrees with the 
expansion at the critical point, where Biq = 1.75 • • •. 

Fig. 4 shows the specific heat variation for the lattices 
with size = 101,201,401,801. Fig. 4a shows the global 
feature in the critical region. The diverging trend as 
t approaches 0 can be clearly seen. In order to show 
this divergence more clearly, we show the part for t > 0 
and t < 0 in Fig. 4b and 4c respectively. The linear 
dependence of the variation on the is shown. As 

the size increases, the linear region grows. Fitting the 
linear parts for N = 800, we get 

(42) 

where C" = 1.8 for t > 0 and C = 1.5 for t < 0. 

As we know, near the critical point it has c ^ 
and || ~ \t\~^ fo'' the two dimensional Ising model [^. 
The finite size scaling near the critical point verifies the 
conclusions given in Eq. (l24l) and (l26ll again. 


IV. SUMMARY AND ACKNOWLEDGEMENT 

We show general properties for the critical systems 
with a single defect. The internal energy variation is 
proportional to the specific heat of the pure system and 
the heat capacity variation is proportional to the tem¬ 
perature derivative of the specific heat. The numerical 
calculation for two dimensional Ising model with a site 
defect verifies these properties. We present a theoretical 
















approach to describe a very small consequence on ther¬ 
modynamic quantities of a defect and their FSS relations 
although the numerical variations are very small 

The effect of a defect is relatively small if the system 
is macroscopic. However if the system size is small, the 
defect’s effect becomes remarkable. Besides, there are 
some advantages in measuring 5U, 6C. 

1. To find the defect’e effect, the relative accuracy 
needs not to be so high. For example, consider a 100 x 
100 lattice (made with cold atom technique). For a non- 
critical system, if one atom is removed, the effect can be 
discovered only if the relative accuracy reaches 1/10000. 
However, if the system is Ising-like and at the critical 
point, the effect can be discovered if the relative accuracy 
of heat capacity reaches about 1/570. This is because the 
total heat capacity is about 0.45 x 10000 x In 100 and the 
variation of the heat capacity is 0.37 x 100 according to 
our numerical result. 

2. Another remarkable feature of 6C is that it changes 
sign as the temperature crosses the critical temperature. 

3. The boundary effects are cancelled in SU because 
[/, Uo have the same boundary effects. It is the same for 
SC. 

Modern nano technology and technology of cold atom 
makes it possible to test these predictions. We are ex¬ 
pecting the experiments on this topic. 

The authors thank Joseph Indekeu, Zehui Deng, 
Adoardo Lauria, Fedorico Galli for useful discussions. 


TABLE IV. The fitted parameters in Eq. (38) for the free 
energy variation with p = 1. It has Smax. < 10“'^®. 

_Ai_ AAj 

Ao 0.1409210310117490418347951D-t01 0.2D-21 
Ai -0.8139392192780112125175834D-t00 0.5D-18 
Aa 0.9067904677646148699557989D-t00 0.5D-15 
As -0.1951243710139618217363658D-t01 0.3D-12 
A 4 0.5044255710564309572237116D-t01 0.9D-10 
As -0.1097266684518414528324830D-t02 0.2D-07 
Ae 0.2276143211351712932688598D-t02 0.3D-05 
At -0.5145819510473408538160129D-t02 0.3D-03 
Ag 0.2111822173858678998575746D-t03 0.2D-01 
Ag -0.8914083610684714743102263D-t03 0.8D-t00 
Aio 0.3184490995226739488807939D-t04 0.2D-t02 
All -0.9963717158837611501676144D-t04 0.4D-t03 
Ai 2 0.3453273003867273296055979D-t05 0.3D-t04 

TABLE V. The fitted parameters in Eq. (38) for the free 
energy variation with p = 3. It has Smax < 

_Ai_ AAj 

Ao 0.1409210310117490418348039D-t01 0.2D-21 
Ai -0.6898025698420019143116894D-t00 0.5D-18 
A 2 0.7246887439974793787485822D4-00 0.5D-15 
As -0.1163435600677590595858805D-t01 0.3D-12 
A4 0.2125184926326692099458813D-t01 0.9D-10 
As -0.4104674530813393591114680D-t01 0.2D-07 
Ae 0.1019038004291499836646201D-t02 0.3D-05 
At -0.2695588248794192948700816D-t02 0.3D-03 
Ag 0.9263301940770832606549970D-t02 0.2D-01 
Ag -0.3477772253050403996019055D-t03 0.9D-t00 
Aio 0.1584171365328676077346061D-t04 0.2D-t02 
All -0.6764747691605833782586485D-t04 0.4D-t03 
Ai 2 0.2185652903261572286215693D4-05 0.3D-t04 


Appendix A: Tables for fitted parameters in Eq. 
(38) 

Table IV-VHI show the fitted parameters in Eq. (38) 
for the free energy variation with p = 1,3, 5, 7,11. 
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TABLE VI. The fitted parameters in Eq. (38) for the free 
energy variation with p — 5. It has Smax < 

_Ai_ AAj 

Ao 0.1409210310117490418348008D+01 0.2D-21 
Ai -0.6895803724829525587022584D+00 0.5D-18 
Aa 0.7253638089447714319588130D+00 0.5D-15 
As -0.1160074030138000567633846D+01 0.3D-12 
A 4 0.2124723799806189308335581D+01 O.lD-09 
As -0.4099892626938691331920487D+01 0.2D-07 
Ae 0.1015677717392225622184707D+02 0.3D-05 
At -0.2690303802284446466680259D+02 0.3D-03 
Ag 0.9253410502505634858180149D+02 0.2D-01 
Ag -0.3474282708150724435713793D+03 0.9D+00 
Aio 0.1580378131891456409396060D+04 0.3D+02 
All -0.6710613069181165806211326D+04 0.4D+03 
Ai 2 0.2149053023432619000704101D+05 0.3D+04 


TABLE VII. The fitted parameters in Eq. (38) for the free 
energy variation with p — 7. It has Smax < 10“^®. 

_Ai_ AAj 

Ao 0.1409210310117490418348016D+01 0.5D-21 
Ai -0.6895799575755920995042567D+00 O.lD-17 
A 2 0.7253669127704436033610930D+00 O.lD-14 
As -0.1160055610993603474765616D+01 0.5D-12 
A 4 0.2124777943031982706300862D+01 0.2D-09 
As -0.4099763470545508107346941D+01 0.3D-07 
Ae 0.1015685122021770300170225D+02 0.5D-05 
At -0.2690305830117266022630290D+02 0.4D-03 
Ag 0.9253334141723530059010477D+02 0.3D-01 
Ag -0.3474247146179242454472337D+03 O.lD+01 
Aio 0.1580295378288385220823113D+04 0.3D+02 
All -0.6709437895591080697242263D+04 0.5D+03 
Ai 2 0.2148323931741508684879155D+05 0.4D+04 


TABLE VIII. The fitted parameters in Eq. (38) for the free 
energy variation with p = 11. It has Smax < 

_Ai_ AAj 

Ao 0.1409210310117490418348135D+01 O.lD-20 
Ai -0.6895799567993295580763121D+00 0.2D-17 
A 2 0.7253669220327549364870848D+00 0.2D-14 
As -0.1160055538390475552667760D+01 0.9D-12 
A 4 0.2124778303511495146132637D+01 0.3D-09 
As -0.4099762126894845643279945D+01 0.5D-07 
Ae 0.1015685523707872014193952D+02 0.7D-05 
At -0.2690310124826811330621585D+02 0.7D-03 
Ag 0.9253638328030327623722488D+02 0.4D-01 
Ag -0.3475505564770179303424906D+03 0.2D+01 
Aio 0.1583637648017210503021256D+04 0.5D+02 
All -0.6760622483924301563712628D+04 0.7D+03 
Ai 2 0.2182556528199001879566179D+05 0.5D+04 


TABLE IX. The fitted parameters in Eq. (40) for the internal 
energy variation with p = 1. It has Smax < 10“^®. 

_ Si _ ABj 

Bio 0.1756000940507335790614587D+01 0.5D-16 
Boo 0.1108789097300574644698700D+01 0.5D-15 
Bn -0.1429278034573753818480900D+01 0.7D-12 
Boi 0.6499233012947408899114641D+00 0.6D-11 
Bi2 0.2173997621082043760717310D+01 0.2D-08 
Bo 2 -0.2578697277246965898741373D+01 O.lD-07 
Bis -0.4426808373674998573819779D+01 O.lD-05 
Bos 0.5993117779650739766451307D+01 0.8D-05 
Bi4 0.1153234993606497236496378D+02 0.4D-03 
Bo 4 -0.1059305245056648116716548D+02 0.2D-02 
Bis -0.2839047458341355573904127D+02 0.4D-01 
Bos 0.2106973174766726014976398D+02 0.2D+00 
Bie 0.2524070175102941639006634D+02 O.lD+01 
Boe 0.9182572691186146023275088D+02 0.5D+01 
Bit -0.7718165760748037716462095D+03 0.2D+02 
Bot 0.1508727615365346140113329D+04 0.3D+02 
Bis -0.1443142955304667511037682D+04 0.4D+02 
Bog -0.1211096405529147735468632D+04 0.2D+02 

TABLE X. The fitted parameters in Eq. (40) for the internal 
energy variation with p = 3. It has Smax < 10“^®. 

_Bi _ ABj 

Bio 0.1756000940507336069227175D+01 0.5D-16 
Boo 0.1442635599708725401635079D+01 0.5D-15 
Bn -0.1211293961408149151708678D+01 0.7D-12 
Boi 0.5043655667225743104254069D+00 0.6D-11 
Bi2 0.1690330956193629515035421D+01 0.2D-08 
Bo 2 -0.1182276963774323003930241D+01 O.lD-07 
Bis -0.2740849022641571930683537D+01 O.lD-05 
Bos 0.1404395814863511902613290D+01 0.8D-05 
Bi4 0.5181119451810307304256970D+01 0.4D-03 
Bo 4 -0.2195324888537194694893548D+01 0.2D-02 
Bis -0.1063018243120927121377483D+02 0.4D-01 
Bos 0.6380322280328616256054328D+01 0.2D+00 
Bie 0.1680697493821303414490378D+02 O.lD+01 
Boe 0.1559103790846971099277808D+02 0.5D+01 
Bit -0.1981485764840919333269028D+03 0.2D+02 
Bot 0.3337347563283782000262181D+03 0.3D+02 
Bis -0.2418971209888612849943305D+03 0.4D+02 
Bog -0.2713954613779996002250673D+03 0.2D+02 

Appendix B: Tables for Eq. (40) 

Table IX-XIII show the fitted parameters in Eq. (40) 
for the internal energy variation with p = 1,3, 5, 7,11. 
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TABLE XL The fitted parameters in Eq. (40) for the internal 
energy variation with p — 5. It has Smax < 

__ ABj 

Bio 0.1756000940507336058650080D+01 0.6D-16 
Boo 0.1443994407066697150440489D+01 0.7D-15 
Bn -0.1210903782636823168187756D+01 0.9D-12 
Boi 0.5091515461841247195439272D+00 0.7D-11 
Bi2 0.1691247267561949045430528D+01 0.2D-08 
Bo2 -0.1161748603765854362514163D+01 0.2D-07 
Bi3 -0.2735654880042078251899909D+01 0.2D-05 
Bo 3 0.1410024855708987543124567D+01 O.lD-04 
Bi4 0.5176567746969383580895496D+01 0.4D-03 
Bo 4 -0.2179708739966420878693034D+01 0.2D-02 
Bis -0.1062524322352045950670683D+02 0.4D-01 
Bos 0.6260914202077016743729399D+01 0.2D+00 
Bi 6 0.1649388105790230667518388D+02 0.2D+01 
Boo 0.1651801479893153943505287D+02 0.5D+01 
Bit -0.2007530400126966705185987D+03 0.2D+02 
Bot 0.3381149739344958924087903D+03 0.3D+02 
Bis -0.2476407744793499744789163D+03 0.4D+02 
Bos 0.2741340607339576860372312D+03 0.2D+02 

TABLE XII. The fitted parameters in Eq. (40) for the inter¬ 
nal energy variation with p — 7. It has Smax < 10~^^. 

_Bi_ ABj 

Bio 0.1756000940507336016832686D-h01 0.3D-15 
Boo 0.1443998353193895611366441D+01 0.3D-14 
Bn -0.1210903054059646131733464D+01 0.3D-11 
Boi 0.5091822951868815096639522D+00 0.3D-10 
Bi2 0.1691252214235173933684408D+01 0.7D-08 
Bo2 -0.1161571472636836367442673D+01 0.5D-07 
Bi3 -0.2735626953648101518756340D+01 0.5D-05 
Bo3 0.1410560164873363940494722D+01 0.3D-04 
Bi 4 0.5176421862069352651115395D+01 O.lD-02 
Bo4 -0.2177402993806319108001286D+01 0.6D-02 
Bis -0.1064664254167213280768083D+02 O.lD+00 
Bos 0.6350090315565562597904054D+01 0.5D+00 
Bi6 0.1571933450600032936027599D+02 0.4D+01 
Boo 0.1886927025564315622330318D+02 O.lD+02 
Bit -0.2092958846293486710359762D+03 0.5D+02 
Bot 0.3525515468089597601811773D+03 0.8D+02 
Bis -0.2646422214482911819703072D+03 0.9D+02 
Bos -0.2846318864144246813607891D+03 0.6D+02 


[1] H. E. Stanley, Rev. Mod. Phys., 71, S358 (1999) 

[2] L. Onsager, Phys. Rev. 65, 117 (1944). 

[3] L. P. Kadanoff, Phys. Rev. Letts., 23, 1430 (1969) . 


Appendix C: Tables for Eq. (41) 

Table XIV-XVIII show the fitted parameters in Eq. 
(41) for the heat capacity variation with p = 1, 3, 5, 7,11. 

TABLE XIII. The fitted parameters in Eq. (40) for the 
internal energy variation with p — 11. It has 5max < 10”^^. 

_ Bi _ ABj 

Bio 0.1756000940507336137644741D-h01 O.lD-14 
Boo 0.1443998363217383095050024D+01 O.lD-13 
Bn -0.1210903052695266140292354D+01 O.lD-10 
Boi 0.5091824170970866970794450D+00 O.lD-09 
Bi2 0.1691252231953949058566169D+01 0.3D-07 
Bo2 -0.1161570544019369750197256D+01 0.2D-06 
Bi3 -0.2735625429719691217440854D+01 0.2D-04 
Bo3 0.1410556584368756569688147D+01 O.lD-03 
Bi4 0.5176729371229483185759416D+01 0.4D-02 
Bo4 -0.2178894119812718899344679D+01 0.2D-01 
Bis -0.1062055808299885422980119D-f02 0.4D+00 
Bos 0.6246745084990453792973751D+01 0.2D+01 
Bi6 0.1654368780944730319389318D+02 O.lD+02 
Boo 0.1646591910591960103095858D-f02 0.4D+02 
Bit -0.2011852559895122908091908D+03 O.lD+03 
Bot 0.3397153318119339570554453D+03 0.2D+03 
Bis -0.2500642654428872245075102D-f03 0.3D+03 
Bos -0.2742309124676317551292031D+03 0.2D+03 

TABLE XIV. The fitted parameters in Eq. (41) for the heat 
capacity variation with p = 1. It has 5max < 10“^^. 

Ci ACi 

C-i 0.3701318932736015169025661D-h00 0.2D-11 
^20 0.5988365596873042415899176D-h00 O.lD-07 
Cio 0.1362181075220199930949431D-01 0.3D-06 
Coo 0.1564365611908562201157505D-1-00 O.lD-05 
^21 -0.9698724835676889360365375D-h00 0.4D-04 
Cii 0.3700976030094435769966437D-h00 0.4D-03 
Coi 0.4738066353822200246665512D-01 0.2D-02 
C 22 0.1369062477761004233681053D-h01 0.4D-02 
C 12 0.1055363190577289429792387D-h01 0.3D-01 
C 02 0.6685542891174151658098970D-1-01 0.7D-01 
C 23 -0.5617895569195019217213644D-h01 0.2D-01 
Cio 0.1968401777398099423864379D-h02 0.3D-01 
Coo -0.4345729961143958935258021D-h01 0.8D-01 


[4] J. Cardy, Scaling and renormalization in statistical 
physics, Cambridge University Press (1996). 

[5] J. Cardy, Conformal Geld theory and statistical mechan¬ 
ics, arXiv 0807.3472 (2008). 
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TABLE XV. The fitted parameters in Eq. (41) for the heat 
capacity variation with p = 3. It has 5max < 10”^®. 

Ci ACi 

C-i 0.4738583363371669290171549D+00 0.2D-10 
C20 0.5988378922929787475482076D+00 O.lD-06 
Cio 0.2412944175798466397089702D+00 0.3D-05 
Coo 0.2240269722942353671837937D+00 O.lD-04 
C21 -0.8253342873928658353818042D+00 0.4D-03 
Cii 0.1680913744352078951227394D+00 0.4D-02 
Coi 0.1173278149931294500887596D+00 0.2D-01 
C22 0.1447584978923885940808479D+01 0.4D-01 
C12 -0.1110998938596696223025486D+01 0.3D+00 
C02 0.1239621370168218395659038D+01 0.7D+00 
C 2 S -0.4878300240223302096521914D+01 0.2D+00 
Cis 0.6192388723449244329513208D+01 0.3D+00 
Cos -0.6071216223290772832895171D+01 0.8D+00 


TABLE XVI. The fitted parameters in Eq. (41) for the heat 
capacity variation with p = 5. It has 5max < 10“^^. 

Ci ACi 

C-i 0.4751207700394731046622277D+00 0.2D-11 
C 20 0.5988378319485564075712419D+00 0.2D-07 
Cio 0.2422223399882923582515648D+00 0.3D-06 
Coo 0.2288063764922809820669021D+00 0.2D-05 
C 21 -0.8249484783990579853710463D+00 0.4D-04 
Cii 0.1693499100245821749277082D+00 0.5D-03 
Coi 0.1401363136723088812062587D+00 0.2D-02 
C 22 0.1437916959377947883645319D+01 0.5D-02 
C 12 -0.1040371814683644478974431D+01 0.3D-01 
C 02 0.1096885820297504552215861D+01 0.8D-01 
C 23 -0.4853343844087464258713497D+01 0.2D-01 
Ci3 0.6288558469821367291050283D+01 0.3D-01 
Co3 -0.5969528111506215066922130D+01 0.9D-01 


TABLE XVII. The fitted parameters in Eq. (41) for the heat 
capacity variation with p = 7. It has 5max < 10~^®. 

Ci ACi 

C-i 0.4751269982626733516613157D+00 0.7D-10 
C20 0.5988375176118030731401894D+00 0.4D-06 
Cio 0.2422310974109595648084283D+00 0.8D-05 
Coo 0.2288248568547370662554431D+00 0.4D-04 
C21 -0.8242508716769945382837724D+00 O.lD-03 
Cii 0.1609297524430012514725075D+00 O.lD-01 
Coi 0.1707615732717604233548380D+00 0.4D-01 
C22 0.1368698698729565558131830D+01 O.lD+00 
C12 -0.6082377549635815123532189D+00 0.6D+00 
C02 -0.4448815393554476152463779D-01 0.2D+00 
C23 -0.4593854244750506845249013D+01 0.4D+00 
Ci3 0.6752451697676445385979327D+01 0.7D+00 
Co3 -0.4709263637062352133372814D+01 0.2D+01 


TABLE XVIII. The fitted parameters in Eq. (41) for the 
heat capacity variation with p — 11. It has Smax < 10~^^. 

Ci ACi 

C-i 0.4751270200380996173992645D+00 0.2D-08 
C 20 0.5988406681372792185253532D+00 O.lD-04 
Cio 0.2421711148698039981400914D+00 0.2D-03 
Coo 0.2291264283210204840237097D+00 O.lD-02 
C 21 -0.8304126144164101519217354D+00 0.3D-01 
Cii 0.2340267004018715595298064D+00 0.3D+00 
Coi -0.8861387725813416786723828D-01 O.lD+01 
C 22 0.1913735962782057699247882D+01 0.2D+01 
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